Dynamical Correlation Functions using the Density Matrix Renormalization Group 
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The density matrix renormalization group (DMRG) method allows for very precise calculations of 
ground state properties in low-dimensional strongly correlated systems. We investigate two methods 
to expand the DMRG to calculations of dynamical properties. In the Lanczos vector method the 
DMRG basis is optimized to represent Lanczos vectors, which are then used to calculate the spec- 
tra. This method is fast and relatively easy to implement, but the accuracy at higher frequencies is 
limited. Alternatively, one can optimize the basis to represent a correction vector for a particular fre- 
quency. The correction vectors can be used to calculate the dynamical correlation functions at these 
frequencies with high accuracy. By separately calculating correction vectors at different frequencies, 
the dynamical correlation functions can be interpolated and pieced together from these results. For 
systems with open boundaries we discuss how to construct operators for specific wavevectors using 
filter functions. 



I. INTRODUCTION 

Since its development, the density matrix renormal- 
ization groupM (DMRG) has been successfully used to 
calculate static properties of ground states and low-lying 
excited states in various low dimensional strongly inter- 
acting systems. Energies can be determined with high- 
est precision, and the calculation of time-independent 
correlation functions is easy and high accuracy can be 
achieved. The calculation of dynamical properties is 
more difficult. 

The Lanczos vector method, also known as the con- 
tinued fractions method, can be used to determine the 
dynamical correlation functions in an exact diagonaliza- 
tion calculation. However, in a DMRG calculation, if 
the basis is optimized only to represent the ground state, 
this will lead to poor results, since the Lanczos vectors 
are not represented correctly in the truncated basis. Hall- 
berg suggested using several of the first Lanczos vectors 
as target states in addition to the ground stateo. Reason- 
able results were obtained for an S=l/2 chain, but the 
true accuracy of the method was not determined, since 
the infinite system DMRG method, rather than the finite 
system method, was used. In this paper we determine the 
accuracy of this method for the more appropriate finite 
system method. 

An alternative approach for generating dynamical 
spectra is the correction vector method, which yields ex- 
act results within the given bjasiai Although it has been 
used in the DMRG contexts, in those calculations the 
basis was not expanded to represent these correction vec- 
tors. Here we apply the correction vector method, target- 
ing a correction vector at a particular frequency. We find 
that the dynamical correlation function at this frequency 
can be calculated directly and very accurately. We show 
how to piece together results from correction vectors at 
different frequencies to obtain the full spectrum. 

Taking antiferromagnetic Hciscnbcrg chains with spin 
1 and spin 1/2 as examples, we discuss the advantages 



and limitations of the two methods. We show how to 
obtain the spectral weight functions, and how to judge 
the quality of the numerical results. 

DMRG calculations are most accurate with open 
boundary conditions, in which case momentum is not 
precisely defined. In this work we show how to construct 
operators corresponding to wavevectors in systems with 
open boundaries using filter functions. 

In Section [l| we discuss the construction of the opera- 



tors for systems with open boundaries. In Section III 



present the Lanczos method, and in Section IV we apply 
it to the antiferromagnetic spin-1 chain. As an example 
where the Lanczos method does not work so well we dis- 
cuss the antiferromagnetic spin-1/2 chain in Section 
In Section VI we present the corr ectio n vector method, 



and we give conclusions in Section VII 



II. CONSTRUCTION OF OPERATORS FOR 
OPEN SYSTEMS 



To calculate a Green's function 

1 



G(q,z) = (0\Al 



z-H 



A q \0) 



(1) 



with z = lo + irj, we have to be able to apply an operator 
A q in our system. The DMRG works in real-space, and 
operators with wavevectors q can be obtained as Fourier 
transforms of on-site operators A n . In infinite systems 
they are given as: 



Aq — ^ ' A n e 



A n — 



2tt 



dqA q e- tx "i 



(2) 



(3) 



where x n is the position of site n, the lattice spacing is 
a = 1. 
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For finite systems with open boundaries, we construct 
operators defined as wavepackets, with finite spatial ex- 
tent and finite uncertainty in the momentum. We con- 
struct the wavepacket by inserting a windowing or filter 
function in Eq.(Q). If only real operators are used, it is 
numerically convenient to construct 



A(q) = ^ sin(qx n )f(x n )A n 

n— — OG 

= ±J dq'(A q ,-A_ q ,)f(q-q>), 



(4) 



and 



Ml) = ^2 cos(qx n )f(x n )A n 

dq'(A q ,+A_ q ,)f(q-q') 



n— — oo 
47T 



(5) 



where f(x n ) = F(jj) is the filter function, 2M is the 
width of the window. We use system with even numbers 
of sites, and x n is offset so that x — is in the center of 
the system. The sites closest to the middle of the system 
are at x — —1/2 and x — 1/2. The operators A(q) are 
reflection symmetric, which allows, if the Hamiltonian is 
also reflection symmetric, using reflected system blocks 
as environment blocks in the DMRG. 

Applying the operator as if the system were periodic, 
is equivalent to using a rectangular window F r (x) as the 
filter function with 2M = L: 



F r (x) 



if -1 < x < 1 
otherwise 



(6) 



This operator is seriously flawed. First, it has sub- 
stantial weight at the edges of the system, where open 
boundary effects are significant. Second, even if we ig- 
nore the edge effects, this window is very broad in the 
wavevector space. The Fourier transform of F r (x) is 



F r (q) 



' Lq 



) (Mi) 



in(g) cos( |) 
cos(q) — 1 



which for small q is 



F r (q) 



2sm(L/2q) 
1 



(7) 



(8) 



The wavevector uncertainty Aq of this operator is of 
order 1, even when L — > oo, whereas q itself ranges from 
to 2ir. Therefore this operator is not useful. 

As is well known from elementary quantum mechan- 
ics, the wavepacket with the minimum product of un- 
certainties AqAx is a Gaussian function. However, it is 
more desirable to use an operator which is equal to at 
the edges of the system. A widely used filter that looks 
similar to a Gaussian in the center, but with compact 
support, is the Parzen filter (Fig. [I]): 

6\x\ 2 + 6\x\ 3 if 0< |af| < 1/2 



F p (x) 



1 



2(1- 



if 1/2 < |ar| < 1 



(9) 
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FIG. 1. Parzen and Gaussian functions in a L = 20 
site system. The halfwidth of the Parzen function is ap- 
proximately 0.18L. The Gaussian has a standard devia- 
tion of a ~ 0.153L, while the standard deviation of the 
Parzen window, because of the faster decay of its tails, is 
only a » 0.112L. 



With 2M = L this filter smoothly goes to zero at the 
boundaries. The Fourier transform is 



F p (q) w 24 



3+cos(qM)-4cos(gM/2) 



(10) 



and the wavevector uncertainty is Aq = 2y3/M. We see 
that Aq varies inversely with the system size if 2M = L. 
For example, in a system with 100 sites a w 0.07, and 
with wavevectors < q < 2ir this is a relative error of 
only Aq — 1.1%. Figure || shows how the Fourier trans- 
form of the Parzen filter smoothly goes to zero, while the 
Fourier transform of a rectangular window has oscilla- 
tions at higher frequencies. 
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FIG. 2. The Fourier transform of the Parzen filter and 
the rectangular window. In real space both niters are 20 sites 
wide. It is obvious that the Parzen filter has a smaller width 
in wavector space. 
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The use of sine (Eq. ||) and cosine (Eq. ||) functions 
for the operators instead of the complex form (Eq. ||) 
requires some special attention. The Green's function is 

= ± Al q )-l^(A q ± A. q )) . (11) 

with "— " for sine and "+" for cosine. Noting that 

( A -9l^ff A -<J> = ( A \j=H A <l) aild A -* = thCTC 

are three different cases for the Green's function. For 
Ml) = En cos(gn)A„ it is: 



\iA\jljjAq) for0<g<7r 
G{q,z) = { {A\^jjA q ) for g = 




(12) 



for q = 7r 



For A(q) = J2 n sin(gn)A„ 



f 5<4i^?A?> for < g < tt 
G(q, z) = { for q = 



(13) 



If only sine or cosine functions are used, the values 
found at q = and q — tt are either zero or twice the 
expected values. By doing separate calculations with sine 
and cosine functions and adding up the results the correct 
values are always obtained. 

In finite systems this effect is broadened due to the fi- 
nite width Aq of the operators in the wavevector space. 
To get system size independent values for the Green's 
functions we require tt 2 = jdqf(q) 2 for the filter func- 
tions. For the Parzen filter this means that an additional 
prefactor of 



1407T 

151A/ 



must be included in the filter. 



III. LANCZOS VECTOR METHOD 

To calculate spectral functions a Lanczos vector pro- 
cedure can be usecB. To do this, the Hamiltonian H is 
projected onto a Krylov subspace spanned by Lanczos 
vectors |/„): 



l/o) 


= A q \0)/{Q\A\A q \Q) 


a n 


= (fn\H\fn) 


b„ 


= \\K)h 


Vn) 


= {H - a n )\f n ) - &„-l|/ n -l 


|/n+l) 


= \r n )/b n 



(14) 



In the Lanczos basis the Hamiltonian is tridiagonal: 

\ 



H 



( ao 
bo 



\ 



b 
a i 
bi 



h 

a-2 



bn-2 a n-l b n -\ 

b n -i a n 



Now the eigenvectors |<& n ) of H are used as an approx- 
imation of the identity 1 w ^ n |$ n )( ( I > n|- Inserting this 
into the Green's function we get: 

G(q,z)«^(0|At|$„)($J-2_|ci> m )($ m |^|0) 

m 

= Y.^\ 7 ^Jj\^n)^n\A q \Q) 2 

n 

_^ {<S>lf{0\A\A q \Q) 

Z £,„ 

Here E n is the eigenvalue of and $™ = (/„|<£„) ■ The 
dynamical correlation function I a u>) is then given by 



Ia{Qi kO = Im lim G(q,u) + irj + E g ) 

TT r)^0+ 

(0\A\A q \0) um ^($0)2 



tt ^ {uj + E g - E n ) 2 + r/ 2 

= (0144,10) J2 S (" + E 9- E n ){K) 2 , (16) 



where B a is the ground state energy. The peaks in the 
correlation function are at u> n = E n — E g . 

The major source of errors in this is the approximation 
of the identity in terms of Lanczos vectors, and the small 
number Nl of these Lanczos vectors that can be used as 
target states in practical calculations. Only those states 
that are used as target states can be relied on as being 
represented correctly. Similar to the Lanczos algorithm 
for the solution of eigensystems, every additional Lanczos 
vector adds another peak, typically at high energy, while 
peaks with the smallest frequencies are determined most 
precisely and converge the fastest. 

To obtain the dynamical correlation function at the 
end of the calculation, we make full use of the available 
Hilbert space. We calculate I(q, u>) at the DMRG step 
where the system block is the same size as the environ- 
ment block, since the truncation errors are smallest in 
that step. To do this, we not only use the Lanczos vectors 
that were target states, but keep calculating new Lanczos 
vectors until orthogonality breaks down and the overlap 
of the new Lanczos vector \f n ) with the first Lanczos 
vector |/o) is bigger than 1%. 

If more than just a few Lanczos vectors are used, the 
question of the weight that is assigned to these target 
states in the density matrix has to be addressed. The 
weight of a Lanczos vector in the spectrum is given by: 



(17) 



Taking this as a measure for the importance of a Lanczos 
vector, we assign 50% of the weight to the ground state, 
and distribute the remaining 50% among the Lanczos 
target states according to their weight uu n . 
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IV. THE SPIN 1 HEISENBERG MODEL 

As an example we look at the antiferromagnetic spin 
1 Heiscnberg model: 



H 



i+1 



(18) 



Each of the spins can be viewed as two spin-1/2 spins 
that pair with the spin-1/2 on the neighboring site in an 
antisymmetric singlet wavefunctionQ. In an open chain 
this effectively leaves unpaired spin-l/2's at the ends. 
To compensate for this, we add real spin-l/2's at the 
ends of the chain. We do not include theses spins in the 
calculation of the spectra and set 2M = L — 2 for the 
Parzen filter. 

This system has a finite correlation length £ w 6.03(1) 
and the Haldane gap Ah = 0.41050(2)0, that separates 
the ground state from the first excitation, a single 
magnon with wavevector tt. DMRG is particularly ac- 
curate in this system even if the number of states kept is 
small. 

To obtain the dynamical correlation function we keep 
the ground state and the three first Lanczos states as tar- 
get states (Nl = 3), and we keep m — 128 states in the 
DMRG basis. To verify that keeping three Lanczos vec- 
tors as target states is enough, we calculate the weight of 
the Lanczos vectors in the final S + (q,uj) for q = tt. Fig- 
ure H shows that the weight of the first vectors is big and 
decays fast, and the first three Lanczos vectors contain 
98.87% of the total weight. 



To further verify the convergence of the DMRG basis, 
we calculate S' + (q,uj) in every step and compare it to 
the final result S + (q, uo) by taking the inner product of 
the two: 



S' + (q) ■ S+(q) 
\\S'+(q)\\\\S+(q)\\ 



(19) 



with the inner product defined as A ■ B = 
f* n du>A(cj)B(uj), and the norm ||A|| = V 'A ■ A. For this 
calculation we use a finite broadening factor r\ = 0.01. 
Figure [| shows the product for every DMRG step. The 
discontinuities occur when the system and environment 
blocks have the same size, and utilizing reflection sym- 
metry, the system block is reflected onto the environment 
block. After two sweeps the DMRG basis is converged, 
there are no further discontinuities, only small oscilla- 
tions which are due to the different truncation effect de- 
pending on the size of the system and the environment 
block. Obviously with three Lanczos target states the 
spectrum is described well enough, and convergence is 
very good. 

In the calculation of the final spectrum, Lanczos vec- 
tors are calculated until the overlap of the new vector 
with the first one exceeds 1%. The inset of Fig. [| shows 
that the overlap (f n \fo) increases exponentially. In the 
given example we use the first 83 Lanczos vectors. It 
can be seen from the overlap that orthogonality breaks 
down after 95 Lanczos vectors have been calculated, ef- 
fectively restarting the procedure. The effect on the re- 
sulting spectrum if more Lanczos vectors are used is very 
small. 



10° 



10" 



10" 



10 E" 



10" 



10" 



• • Spin 1/2 

= Spin 1 



• . • 



20 



40 

n 



60 



80 



FIG. 3. The weight w n of the Lanczos vectors in the 
spectrum. L — 320, Nl = 3 Lanczos target states, m = 128 
states for Spin 1 and m = 256 states for Spin 1/2, q = n. 
Weight for the first three vectors: X/n-o Wn ~ 0-2799 for 
Spin 1/2, and YZ= Wn = 0.9887 for Spin 1. 
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FIG. 4. Deviation from unity of the overlap of the spec- 
tral weight calculated in every step with the spectral weight 
calculated in the last step. Spectral weights at q — n in a 320 
site long spin-1 chain, with a broadening factor rj = 0.01 and 
m = 128 states kept. Iterations are counted from the first step 
after the build-up phase. The dashed lines indicate the iter- 
ations at which the system and environment block have the 
same size. The inset shows the overlap of the first Lanczos 
vector |/o) with \f n ) in the final calculation of the correlation 
function. 
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Figure |5| shows S + (q, ui) for q = ir. Most of the weight 
is in one single peak, and with increasing system size this 
peak moves towards the Haldane gap. To be sure that 
keeping to — 128 states in the DMRG basis is enough, 
we compare results with those from calculations with 
to = 64 and m = 256 states. Table | shows the trun- 
cation error depending on the system size and number of 
states. The truncation errors are relatively small even 
with m = 64 states, and the difference between those 
with to = 128 and to = 256 states are small. 

TABLE I. The truncation error P(m) as a function of the 
length L of the spin-1 chain and the number m of states kept. 
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FIG. 5. The single magnon excitation at q — n in spin-1 
chains of different length versus the frequency minus the Hal- 
dane gap Ah- Calculations with m = 128 states kept and 
broadening factor rj = 0.001 



Table 111 and Fig. |6j(b) show the position of the first 
peak loq. We expect this value to converge against the 
Haldane gap Ah = 0.41050(2) for big systems. For the 
L = 320 site chains our result is only 0.7%(128 states) 
above this value, and Fig. ||(b) indicates that the dis- 
tance between the Haldane gap and the peak vanishes 
for big systems. 

By determining the frequency at which the first peak 
is for different wavevectors, we can calculate the sin- 
gle magnon dispersion relation. Figure [t] shows the 
dispersion relation determined with the Laaczos vector 
method, as well as quantum Monte CarloB and exact 
diagonalizationO. The different sets of data are in good 
agreement, except for small q, where the quantum Monte 
Carlo results are smaller than ours. 

In summary, we have shown that the Lanczos vector 
method works very well for the dynamical spectrum of 
the antiferromagnetic spin-1 Heisenberg model. 
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FIG. 6. The upper plot (a) shows the weight of the first 
peak Wo versus the system size L. The lower plot (b) shows 
the difference of the frequency of first peak u>o and the Hal- 
dane gap Ah ■ The difference tends to zero as the system size 
increases. 



Table || and Fig.|(a) show the weight in the first peak. TABLE III. The position of the first peak lo as a function 

The weight grows with the system size, but seems to go f the length L of the spin-1 chain and number m of states 
to W = 97.6% ± 0.1% instead of 100%. This value is kept, 
important for estimating how good a single mode ap- 
proximation for this excitation is. 



TABLE II. The weight Wo of the first peak as a function 
of the length L of the spin-1 chain and number m of states 
kept. 
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ujo(m — 64) 


Lj {m = 128) 


u>o(m = 256) 
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0.4174 


320 
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L 


Wo(m = 64) 


W (m = 128) 


W (m = 256) 


40 


0.9346 


0.9288 


0.8415 


80 


0.9675 


0.9617 


0.9540 


160 


0.9753 


0.9742 


0.9737 


320 


0.9755 


0.9755 


0.9756 



5 



<4 

3° 



L=1 8 periodic, ED 
L=32 periodic, QMC 
L=80 open, DMRG 
i L=160 open, DMRG 
L=320 open, DMRG 



•fa 



0.0 



0.2 



0.4 



0.6 



0.8 



1.0 



q/n 



FIG. 7. The single magnon line of the Spin 1 Heisenberg 
antiferromagnet. DMRG results with m = 128 states per 
block, and three La»czos vectors as target states Quantum 
Monte Carlo(QMC)El and exact diagonalizationla(ED) data 
are from the work of Takahashi. 



V. THE SPIN-1/2 HEISENBERG MODEL 

As we have already pointed out, the spin-1 chain is a 
relatively easy case. Since there is an excitation gap and 
the correlation length is short, the ground state proper- 
ties can be calculated with high accuracy with a relatively 
small DMRG basis. Following the example of HallbergS, 
we now investigate the antiferromagnetic spin-1/2 chain. 
This is a more difficult case, since it has no excitation gap 
and a diverging correlation length, requiring a slightly 
bigger DMRG basis. More importantly, instead of a sin- 
gle peak with most of the weight in it, for S + (uJjm) an 
excitation band is found. It has a lower boundarytd 



u\ = | I sin(g) | , 



(20) 



and upper boundaryQ 

w» = tt | sin( g /2) | . (21) 

For this band the following structure was proposed^: 
A 



S zz (q,co) 



UJ 2 — UJ l q ' 



=e(u,-^)e«-a,). (22) 



In a finite chain there can be no continuous energy 
band. Instead, we expect separate peaks in the same 
region, and that the number of peaks increases as the 
system size is increased. 

This poses a problem to the Lanczos vector method. 
Figure || shows that the weight of the Lanczos vectors in 
the spectral weight function decreases only very slowly. 
This means that a lot of Lanczos vectors are needed to 
describe the correlation function well. Figure @ shows 



the spectral weight at q — n determined with differ- 
ent numbers of Lanczos vectors. For reference we also 
show results using the much more accurate correction 
vector method, discussed in the next section. The spec- 
tral weight function strongly depends on the number of 
Lanczos target states, and even for 16 and 32 Lanczos 
vector target states there is no sign of convergence. Big- 
ger numbers of target states Nl would result in very long 
calculation times, and would also require that we keep 
substantially more states per block. 
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FIG. 8. The spectral weight function at q = n in a 160 
site spin-1/2 chain with r\ = 0.1 for calculations with different 
numbers of Lanczos target state Nl and with the correction 
vector method. In all calculations m = 256 states were kept in 
the DMRG basis. The correction vector method works well, 
while the Lanczos method is far from convergence. 



While the high energy part of the spectral weight func- 
tion is not accessible with the Lanczos vector method, 
the position of the first peak is only weakly dependent 
on the number of Lanczos vector target states Nl . Table 
[V shows the position of the first peak wo in a L = 160 
site chain at q = n calculated with m = 256 states kept. 
The dependence of the position of the first peak on the 
number of Lanczos vector target states is small. The rea- 
son why it is shifted to higher values if more target states 
are used are the increasing truncation errors arising from 
targeting more states with fixed m — 256. 



TABLE IV. The position of the first peak coo depending on 
the number of Lanczos vector target states Nl in the spectral 
weight function of an antiferromagnetic L = 160 site spin-1/2 
chain at q — n. m = 256 states were kept in the DMRG basis. 
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FIG. 9. The lower bound of the excitation band in the 
Spin 1/2 Heisenberg model. The results for the 28 site peri- 
odic system are from Hallbergp. 

Keeping m — 256 states in the DMRG basis, and 
targeting four Lanczos vectors, we have determined the 
lower edge of the energy band for different q. Figure ^ 
shows the dispersion relation. For the longer systems it 
compares increasingly well with the analytic result (Eq. 

E§). 

VI. CORRECTION VECTORS 

In the previous section we found that the Lanczos vec- 
tor method works very well if only low-energy proper- 
ties of the dynamical correlation function are of interest, 
but that it is unable to reproduce higher energy proper- 
ties like the shape of the excitation band in the spin-1/2 
model. Instead of using the Lanczos vector method, the 
spectrum can be calculated directly for a given z = Lo + irj 
by using a correction vector. To do this, the following 
states must be included as target states: 

|0) the ground state 

\A q ) = A q |0) the first Lanczos vector 

\ x ( z )) = 7=77 l^<?) the correction vector 

Since a finite broadening factor 77 is used, the correc- 
tion vector becomes complex. To avoid the use of com- 
plex numbers, we split the correction vector into real and 
imaginary part, both used as target states: 

\x(z)) = \x r (z))+i\x i (z)). (23) 

The equation for the correction vector is split into real 
and imaginary parts \x r (z)) and |a; l (z)), and the imagi- 
nary part is found by solving 

{{H-uof+ v 2 ) \x i (z))=-r ] \A q ) (24) 

using the conjugate gradient method. Note that this 
equation system gets more singular as lo gets closer to 



an eigenenergy of the Hamiltonian, and as the broaden- 
ing factor rj gets smaller. For large r\ the Hamiltonian 
is close to a diagonal matrix, and the conjugate gradi- 
ent method converges much faster than for small r\. This 
means that a large rj results in short calculation times, 
but it also limits the resolution of the spectrum. The 
convergence is also slowed down because energy gaps in 
H — lo are squared in Eq. (|24|) , and the convergence rate 
of the conjugate gradient method depends on the gap 
between the lowest and the next lowest eigenvector. 

The real part of the correction vector is calculated di- 
rectly: 

\ X r{z)) = hu-H)\x i {z)) ] . (25) 
V 

Using the correction vector, the Green's function can 
be calculated directly: 

G(q,z) = (A g \x(z)) (26) 

Taking these states (|0),|A g ) and \x(z))) as target states 
and optimizing the DMRG basis to represent them al- 
lows for a very precise calculation of the Green's func- 
tion for a given frequency to and broadening factor rj. 
Unfortunately, the correction vector has to be calculated 
separately for every to. If the correction vector does not 
change very rapidly with lo, the DMRG basis that is op- 
timized to represent the correction vector at a certain 
lo, should also be able to represent correction vectors for 
nearby frequencies. 

Although the correction vectors are needed as target 
states to determine the DMRG basis, it is more efficient 
to use the Lanczos vector method to determine the com- 
plete spectra within that basis. Using the correction vec- 
tor, but no Lanczos vectors except the first one, as a tar- 
get state, DMRG sweeps are performed until the basis is 
converged (typically two or three sweeps). The dynami- 
cal correlation function is then calculated in the same way 
as in the Lanczos vector method: when the left and right 
block have the same size Lanczos vectors are calculated 
until orthogonality is lost. Since this method is almost 
exact in the given basis, in principle it should yield the 
same result as using the correction vector method in the 
same basis, and it is much faster. Of course, while the 
spectrum is produced for all to, it is only accurate near 
the to used to produce the correction vector. 

To determine the range of lo in which the truncated 
basis is good enough to calculate the spectral function, 
we compare calculations with different frequencies for the 
correction vectors. Figure shows the spectral weight 
function near the upper edge of the excitation band in the 
antiferromagnetic spin-1/2 chain. This is an especially 
difficult region, because there are a lot of peaks at lower 
energies. We did several calculations with correction vec- 
tors for different to, calculating the spectral weight in the 
region around the targeted frequencies. With m = 128 
states kept, the overall shape reproduces the edge of the 
band accurately. Deviations between the pieces of the 
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spectrum calculated with different correction vectors can 
be seen. The small mismatch between the very accurate 
values calculated directly from the correction vectors, 
and the corresponding part obtained with the Lanczos 
vector method are due to numerical errors in the Lanc- 
zos procedure. By keeping more states this result can be 
improved. Figure [ll] shows the same part of the spec- 
trum with m = 256 states kept. Now the different parts 
match very well. Further improvement would be possible 
by keeping more states or reducing the distance between 
the frequencies of the correction vectors. 
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FIG. 10. The spectral weight in a 80 site spin-1/2 chain 
at q — n with one correction vector as a target state and 
128 states in the DMRG basis and r\ — 0.1. The crosses show 
the results calculated directly from the correction vectors, the 
lines show the parts calculated with the Lanczos procedure in 
the basis optimized for the correction vectors. For example, 
the dashed line from cj = 2.8totj = 3.0is calculated in the 
basis that is optimized for the correction vector with lo = 2.9. 
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FIG. 11. The spectral weight in a 80 site spin-1/2 chain 
at q = 7r with one correction vector as a target state and 256 
states in the DMRG basis and r\ — 0.1. 



This result can be improved even further. Instead of 
using just one correction vector as a target state, we try 



using two correction vectors at the same time. The spec- 
tral weight can then be interpolated for frequencies be- 
tween these two correction vectors. With a broadening 
factor oirj — 0.1, a distance of Alo = 0.2 between the two 
correction vectors seems appropriate. In Fig. [l^ and [l3| 
the results from calculations with different frequencies 
for the correction vectors are plotted. With m = 128 
states the parts for 3.0 < u> < 3.3 match perfectly, and 
for 2.8 < u> < 3.0 they still match better than with only 
one correction vector as a target state (Fig. |l0|). With 
77i = 256 states (Fig. |l3|) all pieces of the spectrum match 
up almost perfectly. This gives us a consistent method 
to verify how good our numerical results are, and we find 
that it is possible to achieve very high accuracy. 
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FIG. 12. The spectral weight function in a 80 site spin-1/2 
chain at q — n, with 128 states in the DMRG basis and 
rj — 0.1. Two correction vectors are used as target states, 
and the plot shows the values calculated directly with these 
correction vectors, and the connecting lines show the interpo- 
lated spectral weight calculated with the Lanczos procedure 
in the basis optimized for the correction vectors. 



Now we use this method to obtain the complete spec- 
tral weight in a 160 site chain. Keeping m = 256 states 
and using two correction vectors as target states, we start 
with correction vectors for w = and lo = 0.2. Af- 
ter two sweeps through the system the DMRG basis is 
converged, and the spectral weight function is calculated 
for < lo < 0.2 using the Lanczos procedure. Then 
we target u> = 0.2 and lo — 0.4, perform two sweeps 
through the system, and calculate the spectral weight for 
0.2 < lo < 0.4. Continuing this procedure, we obtain the 
function shown in Fig. |[ In contrast to the results from 
the Lanczos vector method, the spectral weight function 
calculated using two correction vectors shows no finite- 
size peaks and reflects the shape expected in the infinite 
system. Eq. ( p2|) predicts 1/lo decay from lo = to 
lo = 7T, where the spectral weight drops to zero. The 
band presented in Fig. |^ shows the correct upper and 
lower bound, but the spectral weight in the band decays 
faster than 1/lo. To verify if this is a finite-size effect, 



8 



we look at chains with different lengths. Figure |IJ shows 
that the spectral weight decays faster than X/uj in all 
system, and the decays do not become slower for longer 
chains. The only hint at hnite-size effects is the different 
size of the spectral weight for different system sizes. 

The small oscillations in the 40 site chain are due to the 
limited number of peaks in the relatively small system, 
those in the 160 site chain are due to the large distance 
of Auj = 0.4 between the correction vectors. They do not 
affect our result, and could be removed by using pairs of 
correction vectors with closer frequencies or by increasing 
the size of the DMRG basis. 
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FIG. 13. The spectral weight function in a 80 site spin-1/2 
chain at q — n, with 256 states in the DMRG basis and 
77 = 0.1. Two correction vectors are used as target states, 
and the plot shows the values calculated directly with these 
correction vectors, and the connecting lines show the interpo- 
lated spectral weight calculated with the Lanczos procedure 
in the basis optimized for the correction vectors. 



In Fig. [13] the upper edge of the excitation band can 
be seen, but due to the broadening factor 77 = 0.1 it is 
difficult to determine the position of the edge. Since the 
spectral weight function is plotted by taking the peaks 
found in the Lanczos procedure and plotting them with 
the given 77, they can as well be plotted with smaller 
broadening factors. Figure |Te| shows plots with 77 = 0.1, 
for which the DMRG basis was optimized, and 77 = 0.05 
and 77 = 0.01. For the smaller rj the different parts of the 
spectral weight function no longer match, but now the 
poles are easily identifiable. The position of the peaks 
can also be determined by directly taking the position of 
the peaks from the Lanczos procedure. Doing this, the 
last peak in the band is found at u> = 3.128, compared to 
uj = n given in (Eq. |2l]). To verify if the small difference 
is a finite-size effect or a numerical error, longer chains 
can be studied, or correction vectors closer to the desired 
region could be used. Using smaller 77 and reducing the 
distance between the correction vectors also gives higher 
resolution. 
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FIG. 14. The spectral weight function at q — n, with 
different system sizes. The spectra were calculated with two 
correction vector target states with distances of Auj = 0.4, a 
broadening factor rj = 0.2 and m = 256 states. 



FIG. 15. The spectral weight function in a 80 site chain 
at q = n. This is the same data as in figure with two 
correction vectors with a distance of Auj = 0.2 between them, 
and 77 = 0.1. If the spectral weight function is plotted with 
a smaller rj, the parts calculated with different frequencies 
for the correction vectors no longer match perfectly, but the 
peaks become easily distinguishable. 



VII. CONCLUSIONS 

In conclusion, we have presented and tested two meth- 
ods to calculate dynamical correlation functions using 
DMRG. We have shown that the Lanczos vector method 
works very well if only the low-energy part of the correla- 
tion function is of interest, or if the bulk of the weight is 
in one single peak. We demonstrated this for the case of 
the antiferromagnetic spin-1 and spin-1/2 chains, where 
we could reproduce the dispersion relation of the lower 
edge of the spectral weight functions correctly. 
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If there is an excitation band, the Lanczos vector 
method is unable to describe the higher energy part of 
the correlation functions. These parts can be determined 
using the correction vector method. This method gives 
very precise results for frequencies at which the correction 
vectors are used as target states. The Lanczos procedure 
can be used in the basis optimized for the correction vec- 
tors to determine the spectrum fast and efficiently not 
only at the frequency of the correction vector, but also 
in the region around that. By comparing the plots from 
calculations with different frequencies for the correction 
vectors it is possible to estimate the range over which the 
spectral function is determined correctly. We find that 
remarkably good spectra can be determined if two cor- 
rection vectors are used as target states, and the spectral 
function is calculated for the frequencies between them. 
In the case of the excitation band in the spectral weight 
function of the antiferromagnetic spin- 1/2 chain, we used 
this method and were able to study a system long enough 
and with sufficient accuracy that no finite-size peaks were 
visible and excellent agreement with theoretical expecta- 
tions was obtained. 

Our results show that using these techniques, obtain- 
ing accurate dynamical spectral functions from DMRG is 
feasible and can be considered a standard DMRG tech- 
nique. The calculation time is longer than for ground 
state properties, but still manageable for single chain sys- 
tems and probably for ladders with a few chains. 
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